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Abstract 

Binding interactions between proteins and other molecules mediate numerous cellu- 
lar processes, including metabolism, signaling, and regulation of gene expression. These 
interactions evolve in response to changes in the protein's chemical or physical envi- 
ronment (such as the addition of an antibiotic), or when genes duplicate and diverge. 
Several recent studies have shown the importance of folding stability in constraining 
protein evolution. Here we investigate how structural coupling between protein folding 
and binding - the fact that most proteins can only bind their targets when folded - 
gives rise to evolutionary coupling between the traits of folding stability and binding 
strength. Using biophysical and evolutionary modeling, we show how these protein 
traits can emerge as evolutionary "spandrels" even if they do not confer an intrinsic 
fitness advantage. In particular, proteins can evolve strong binding interactions that 
have no functional role but merely serve to stabilize the protein if misfolding is dele- 
terious. Furthermore, such proteins may have divergent fates, evolving to bind or not 
bind their targets depending on random mutation events. These observations may ex- 
plain the abundance of apparently nonfunctional interactions among proteins observed 
in high-throughput assays. In contrast, for proteins with both functional binding and 
deleterious misfolding, evolution may be highly predictable at the level of biophysical 
traits: adaptive paths are tightly constrained to first gain extra folding stability and 
then partially lose it as the new binding function is developed. These findings have 
important consequences for our understanding of fundamental evolutionary principles 
of both natural and engineered proteins. 

Proteins carry out a diverse array of chemical and mechanical functions in the cell, 
ranging from metabolism to signaling [1]. Thus proteins serve as central targets for natural 
selection in wild populations, as well as a key toolbox for bioengineering novel molecules 
with medical and industrial applications [2, 3]. Most proteins must fold into their native 
state, a unique three-dimensional conformation, in order to perform their function, which 
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typically involves binding a target molecule such as DNA, RNA, another protein, or a small 
ligand [1]. Misfolded proteins may also form toxic aggregates and divert valuable protein 
synthesis and quality control resources [4-7] . It is therefore imperative that the folded state 
be stable against the thermal fluctuations present at physiological temperatures. However, 
biophysical experiments and computational studies reveal that most random mutations in 
proteins destabilize the folded state [8, 9], including mutations that improve function [9-11]. 
As a result many natural proteins tend to be only marginally stable, mutationally teetering 
at the brink of substantial unfolding [12, 13]. With proteins in such a precarious evolutionary 
position, how can they evolve new functions while maintaining sufficient folding stability? 

Directed evolution experiments have offered a window into the dynamics of this pro- 
cess [2, 3], indicating the importance of compensatory mutations, limited epistasis, and 
mutational robustness. Theoretical efforts to describe protein evolution in biophysical terms 
have focused on evolvability [14], global properties of protein interaction networks [15, 16], 
and reproducing observed distributions of protein stabilities and evolutionary rates [13, 17- 
19]. However, a subtle but key property of proteins has not been explored in this context: 
structural coupling of folding and binding (the fact that folding is required for function) 
implies evolutionary coupling of folding stability and binding strength. Thus selection act- 
ing directly on only one of these traits may produce apparent, indirect selection for the 
other. The importance of this effect was popularized by Gould and Lewontin in their in- 
fluential paper on evolutionary "spandrels" [20] , defined as traits that evolve as byproducts 
in the absence of direct selection. Since then the importance of coupling between traits 
has been explored in many areas of evolutionary biology [21], including various molecular 
examples [12, 22, 23]. 

How do coupled traits affect protein evolution? We consider a simple model that describes 
evolution of a new binding interaction in the context of a directed evolution experiment [3] , 
as a result of gene duplication and divergence [24], or in response to a change in the pro- 
tein's chemical or physical environment, including availability and concentrations of various 
ligands [25, 26] as well as temperature [27, 28]. We postulate a fitness landscape as a func- 
tion of two biophysical traits: stability and the free energy of binding a target molecule. 
We then use an exact numerical algorithm [29, 30] to quantitatively characterize adapta- 
tion on this fitness landscape, addressing key evolutionary questions of epistasis [31, 32], 
predictability [25, 33, 34], and the tempo of adaptation [17, 35]. 

Results 

Model of protein energetics. We consider a protein with two-state folding kinetics [1]. 
In the folded state, the protein has an interface that binds a target molecule. Because the 
protein can bind only when it is folded, the binding and folding processes are structurally 
coupled. Under the thermodynamic equilibrium assumption (valid when protein folding and 
binding are faster than typical cellular processes), the probabilities of the three structural 
states - folded and bound (f>f,b), folded and unbound (pf, u b), and unfolded and unbound 
(Puf,ub) _ are given by their Boltzmann weights: 
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Here /3 is the inverse temperature, Ef is the free energy of folding (also known as AG), and 
Eb = E f b — /x, where E' b is the binding free energy and /i is the chemical potential of the 
target molecule. For simplicity, we will refer to E\> as the binding energy. Note that Ef < 0 
for intrinsically-stable proteins and < 0 for favorable binding interactions. The partition 
function is Z = e -P& f +E b ) + e -pE f + 1 

The folding and binding energies depend on the protein's genotype (amino acid sequence) 
a. We assume that adaptation only affects "hot spot" residues at the binding interface [36, 
37]; the rest of the protein does not change on relevant time scales because it is assumed 
to be already optimized for folding. If positions away from the binding interface can accept 
stabilizing mutations (and are not functionally constrained), they may be explicitly included 
into the model as "folding hotspots." In the present study we focus on L binding hotspot 
residues which, to a first approximation, make additive contributions to the total folding 
and binding free energies [38] (see SI Methods for the discussion of non-additive effects): 

L L 

E f (a) = Ef + °% E b (a) = Ef" + £ e b (i, a% (2) 

1=1 1=1 

where e/(i, a 1 ) and e&(i, a 1 ) capture the energetic contributions of amino acid a 1 at position i. 
The reference energy Ej ei is the fixed contribution to the folding energy from all other residues 
in the protein. Furthermore, by construction it is also the total folding energy of a reference 
sequence a re f (see Methods), so that each €f(i,a l ) can be interpreted as the change in total 
folding free energy Ef (AAG value) resulting from a single-point mutation of cr re f. The 
parameter E™ in is the minimum binding energy among all genotypes (see Methods). Amino 
acid energies e/(i, a 1 ) and e&(i, a 1 ) are randomly sampled from distributions constructed using 
available AAG data and other biophysical considerations (see Methods); the exact shape of 
these distributions is unimportant for large enough L due to the central limit theorem. 

Fitness landscape. We construct a simple fitness landscape based on the molecular 
traits Ef and E b . Without loss of generality, we assume that the protein contributes fitness 
1 to the organism if it is always folded and bound. Let / u b, f u f £ [0, 1] be the multiplicative 
fitness penalties for being unbound and unfolded, respectively: the fitness is / u b if the protein 
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is unbound but folded, and / u b/uf if the protein is both unbound and unfolded. Then the 
fitness of the protein averaged over all three possible structural states in Eq. 1 is given by 

T{E h E b ) = p f5b + iubfVb + /ub/ufPuf,ut>. (3) 

This fitness landscape is divided into three nearly-flat plateaus corresponding to the three 
protein states of Eq. 1, separated by steep thresholds corresponding to the folding and binding 
transitions (Fig. 1A). The heights of the plateaus are determined by the values of / u b and 
/ u f, leading to three qualitative regimes of the global landscape structure (Fig. 1B-D). 

In the first case (Fig. IB), a protein that is perfectly folded but unbound has no fitness 
advantage over an unbound and unfolded protein: / u b = /ub/uf- Thus selection acts directly 
only on the binding trait. This regime requires that either / u b = 0 (binding is essential, 
e.g., in the context of conferring antibiotic resistance to the cell [25]) or / u f = 1 (misfolded 
proteins are not toxic). The latter case also includes directed evolution experiments where 
only function is artificially selected for in vitro. In contrast, when / u b = 1 and 0 < / u f < 1 
(Fig. 1C), a perfectly folded and bound protein has no fitness advantage over a folded but 
unbound protein, and thus this case entails direct selection only for folding. These proteins 
are harmful to the cell in the misfolded state (e.g., due to aggregation or significant costs 
of degrading unfolded proteins [4-7]), while binding provides no intrinsic fitness advantage 
(the protein may have other, functional binding interfaces). Finally, it is also possible to 
have distinct selection pressures on both binding and folding. This occurs when 0 < / u b < 1 
and 0 < /uf < 1 (Fig. ID). 

It is straightforward to generalize our three-state model to proteins with additional struc- 
tural states (other local minima on the folding energy landscape, other binding modes) 
and allow for simultaneous adaptation at multiple binding interfaces. Furthermore, the fit- 
ness landscape in Eq. 3 can be made an arbitrary nonlinear function of state probabilities. 
However, these more complex scenarios would still share the essential features of our ba- 
sic model: coupling between folding and binding traits and sharp fitness thresholds between 
bound/unbound and folded/unfolded states. Thus our qualitative conclusions do not depend 
on the specific model in Eq. 3. 

Epistasis and local maxima. For protein sequences of length L with an alphabet of 
size fc, each of the k L possible genotypes is projected onto the two-dimensional trait space 
of Ef and E\> (Eq. 2) and connected to L(k — 1) immediate mutational neighbors, forming a 
network of states that the population must traverse (a simple example is shown in Fig. IE). 
Adaptive dynamics are determined by the interplay between the structure of the fitness 
landscape and the distribution of genotypes in trait space. 

This interplay gives rise to the possibility of epistasis and multiple local fitness maxima. 
Our model is non-epistatic in energy space (Eq. 2). When the fitness contours are straight 
parallel lines, there can be no sign epistasis on the fitness landscape (Fig. IF). Magnitude 
epistasis, on the other hand, is widespread due to the nonlinear dependence of fitness on 
folding and binding energies. Curved fitness contours, which occur near folding or binding 
thresholds in our model (Fig. 1B-D), can produce sign epistasis in fitness, giving rise to 
multiple local fitness maxima in the genotype space (Fig. IE). 
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Evolutionary dynamics. We assume that a population encoding the protein of interest 
evolves in the monomorphic limit: LNulogN <C 1, where L is the number of residues, TV 
is an effective population size, and u is the per-residue probability of mutation per genera- 
tion [39] (see SI Methods). In this limit, the entire population has the same genotype at any 
given time, and the rate of substitution from the current genotype to one of its mutational 
neighbors is given by Eq. 1 in SI Methods. We use the strong-selection limit of the substi- 
tution rate (Eq. 2 in SI Methods), in which the effective population size enters only as an 
overall time scale. In this regime, deleterious mutations never fix and adaptive paths have 
a finite number of steps, terminating at a global or local fitness maximum. For compact 
genomic units such as proteins, the monomorphic condition is generally met in multicellular 
species, although it may be violated in some unicellular eukaryotes and prokaryotes [40]. Se- 
quential fixation of single mutants is also a typical mode of adaptation in directed evolution 
experiments [3]. For simplicity, we neglect more complex mutational moves such as indels 
and recombination. 

Far from the binding and folding thresholds the fitness landscape becomes flat (Fig. 1A) 
and the strong-selection assumption may be violated. To establish the limits of validity 
for our model, we calculate average selection coefficients of accessible substitutions (defined 
as s = ^final/initial ~~ 1? where initial and iinai are the initial and final fitness values of a 
substitution), both throughout the landscape and at the local maxima (Fig. SI). We observe 
that for typical values of the effective population size TV e (10 4 , 10 7 ) [40, 41], the selection 
strengths in the model justify our strong-selection approximation for realistic choices of 
energy parameters. 

Quantitative description of adaptation. Although our model is valid for any adap- 
tive process, for concreteness we focus on a specific but widely-applicable scenario. A popu- 
lation begins as perfectly adapted to binding an original target molecule characterized by an 
energy matrix with minimum binding energy E™ n (defining a fitness landscape i). The 
population is then subjected to a selection pressure which favors binding a new target, with 
energy matrix e^ 2 and minimum binding energy E™ m (fitness landscape T2). The adaptive 
paths are first-passage paths leading from the global maximum on T\ to a local or global 
maximum on J^, with fitness increasing monotonically along each path. 

Each adaptive path tp with probability U[<p] is a sequence of genotypes connecting initial 
and final states. Using an exact numerical algorithm (SI Methods) [29, 30], we determine the 
path-length distribution p(^), which gives the probability of taking an adaptive path with 
t amino acid substitutions, and the mean adaptation time i. We also introduce /S pat h, the 
entropy of the adaptive paths: 



The path entropy is maximized when evolution is neutral, resulting in all paths of a given 
length being accessible and equally likely: S pat h = l\og L{k — 1) [30], where ^is the average 
path length. 

We also consider the path density ^{cr) ) which gives the total probability of reaching a 




(4) 
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state a at any point along a path. When a is a final state (a local fitness maximum on 
the path density is equivalent to the commitment probability. We calculate the entropy S com 
of the commitment probabilities as 

-Scorn = " ^)l0g^%)- (5) 

final states a 

Direct selection for binding only. We first focus on the / ub = / u f/ u b case i n Eq. 3. 
The geometry of the fitness contours is invariant under overall shifts in the binding energy 

(Fig. IB); equivalently, the direction (but not the magnitude) of the selection force 
(Vlog J 7 / 1 V log J 7 !) does not depend on E^. Thus without loss of generality, we set E™ n = 
E™ n in this section. The contours of constant fitness are parallel to the Ef axis when Ef is 
low, indicating that, as expected, selection acts only on binding when proteins are sufficiently 
stable. 

However, for marginally stable proteins [12, 13, 42], the fitness contours begin to curve 
downward, indicating apparent, indirect selection for folding, even though selection acts 
directly only on the binding trait. Thus, adaptation will produce a trait (more stability) 
that is neutral at the level of the fitness function simply because it is coupled with another 
trait (binding) that is under selection. Folding stability can therefore be considered an 
evolutionary spandrel [20]. Proteins may even be intrinsically unstable (Ef > 0) and only 
fold when bound (Ef + Ef, < 0), which we refer to as binding-mediated stability [43]. In this 
regime, the fitness contours approach diagonal lines: selection effectively acts to improve 
both binding and folding equally (Fig. IB). 

An example realization of evolutionary dynamics in the marginally stable regime is shown 
in Fig. 2A,B (see Fig. S2 for stable and intrinsically unstable examples, and Fig. S3 for 
averaged distributions of initial, intermediate, and final states). There is typically just one 
or two fitness maxima; all maxima are usually accessible (Fig. 2C). For stable proteins, the 
global maximum almost always coincides with the best-binding genotype and is usually as 
far as a randomly-chosen genotype from the best-folding genotype (Fig. 2D; two random 
sequences are separated by 1 — 1/k = 0.8 for k — 5). However, as Ef becomes greater, 
the average distance between the maxima and the best-binding genotype increases, while 
the average distance between the maxima and the best-folding genotype decreases, until 
they meet halfway for intrinsically unstable proteins, where effective selection for binding 
and folding is equally strong (Fig. 2D). In general the maxima lie on or near the Pareto 
front [44], defined here as the set of genotypes such that either Ef or cannot be decreased 
further without increasing the other (the global maximum is always on the front, while local 
maxima may not be) (Fig. 2A, Fig. S2). 

As Ef increases, the average distance between initial and final states for adaptation 
decreases. As a result the average path length (number of substitutions) decreases as well, 
although the variance of path lengths is relatively constant over all energies (Fig. 2E). The 
path entropy per substitution Sp^/l also decreases with Ef, reflecting greater constraints 
on adaptive paths (note that S pat h/^ = log L(k — 1) « 3.2 for neutral evolution). Finally, 
Scom ~ 0.31 in the marginally stable regime (Fig. 2F). Since the average number of maxima 
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is « 1.9 in this regime (Fig. 2C), the maximum value of S com is log 1.9 « 0.64, indicating 
that not all maxima are equally accessible. 

Direct selection for folding only. In this regime, / ub = 1 and 0 < / u f < 1 in Eq. 3. 
Similar to the previous case, the geometry of the fitness contours and thus most landscape 
properties are now independent of Ef (Fig. 1C); equivalently, normalized selection force 
Vlog J 7 / 1 V log J 7 ! does not depend on Ef. 

When the nonfunctional binding is weak, the fitness contours are parallel to the E^ axis, 
indicating that selection acts only on folding (Fig. 1C). However, with increasing binding 
strength the fitness contours curve such that the effective selection force attempts to improve 
both binding and folding equally. Thus binding emerges as an evolutionary spandrel in this 
case. The weak-binding regime yields a single fitness maximum due to the lack of sign 
epistasis; this maximum predominantly coincides with the best-folding genotype (Fig. 3A). 
However, once the binding interaction becomes stronger, there is an increased likelihood of 
multiple local maxima, located between the best-folding and best-binding genotypes. 

Depending on the abundance of the old and new ligands in the cell and their binding 
properties, several adaptive scenarios may take place. First, the best-binding strengths E™ n 
and E™ n of the old and new targets may be similar in magnitude. If both are weak, initial 
and final states are likely to be the best-folding genotype or close to it (Fig. 3A); in this 
case, there is a high probability that no adaptation will occur (Fig. 3B). When E™ n and 
E™ n are both low, adaptation usually occurs to accommodate the binding specificity of the 
new ligand (Fig. 3B, Fig. S4A). Surprisingly, we see that proteins frequently evolve stronger 
binding at the expense of folding (bottom panel of Fig. S4A). This happens due to the 
constraints of the genotype-phenotype map: not enough genotypes are available to optimize 
both traits simultaneously. 

It is also possible to gain or lose binding affinity at the nonfunctional interface through 
adaptation. In the first case, the new target has stronger binding than the old one (Ej£ m < 
E™ n ). Thus the initial state is the best-folding genotype or close to it, and the protein adapts 
toward a genotype with intermediate folding and binding (Fig. S4B). As before, adaptation is 
tightly constrained by the genotype-phenotype map, sacrificing the trait (folding stability) 
under direct selection in order to affect the spandrel (nonfunctional binding interaction). 
Effectively, the protein switches from being "self-reliant" to needing a binding partner. In 
the second case (E™ n < E^ n ) ) the dynamics is opposite: the protein loses its nonfunctional 
binding interface and becomes self-reliant (Fig. S4C). Thus proteins may acquire or lose 
binding interfaces depending on the availability of ligands that can participate in binding- 
mediated stability. If the protein's stability becomes suboptimal due to an environmental 
change, its stability may be restored not only through stabilizing mutations, but also by 
developing a novel binding interface. 

Divergent evolutionary fates. In the region where the fitness contours in Fig. 1C are 
curved, it is possible to have two or more local maxima accessible to adaptation, with at 
least one having negative E\> (strong binding) and at least one having positive Ef, (negligible 
binding) (see Fig. 3C,D for an example landscape). The selection streamlines are divergent 
in this regime (Fig. 1C). Thus a protein has two fates available to it: one in which it evolves 
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to bind the target and another in which it does not. The eventual fate of the protein is 
determined by random mutation events. Indeed, the distribution of final states is strongly 
bimodal (Fig. 3E), yielding a sizable probability of divergent fates across a range of binding 
energies (Fig. 3F). 

Simultaneous selection for binding and folding. Finally we consider a general 
case in which 0 < / u b < 1 and 0 < / u f < 1 in Eq. 3 (Fig. ID). The fitness landscape 
is divided into two regions by a straight diagonal contour with fitness / u b and slope — 1. 
Below this contour, the landscape is qualitatively similar to the case of selection for binding 
only (Fig. IB), while above the contour the landscape resembles that of the folding-only 
selection scenario (Fig. 1C). Thus evolutionary dynamics for proteins with favorable binding 
and folding energies will largely resemble the case of selection for binding only. However, 
a qualitatively different behavior will be observed if the distribution of genotypes straddles 
the diagonal contour (Fig. 4). This will occur when initial folding stability is marginal and 
initial binding is unfavorable. In this case, selection streamlines around the diagonal contour 
(Fig. ID) and the genotype-phenotype map tightly constrain the adaptive paths to gain 
extra folding stability first, and then lose it as the binding function is improved. 

Tempo and rhythm of adaptation. The strength of selection is the primary deter- 
minant of the average adaptation time t. If the selection coefficient s is small (but Ns > 1), 
the substitution rate W^'la) in SI Methods Eq. 1 is proportional to s. Thus, as selection 
becomes exponentially weaker for lower energies (Fig. SI), adaptation becomes exponen- 
tially slower. The distribution of the total adaptation time over an adaptive path is highly 
nonuniform. For example, in the case of selection for binding only and a marginally stable 
protein, the adaptation time is concentrated at the end of the path, one mutation away 
from the final state (Fig. S5A,B). Substitutions at the beginning of the path occur quickly 
because there are many possible beneficial substitutions and because selection is strong; in 
contrast, at the end of the path adaptation slows down dramatically as beneficial mutations 
are depleted and selection strength weakens. This behavior is observed in most of the other 
model regimes as well. 

The exception to this pattern occurs in the case of selection for both binding and folding 
in marginally-stable and marginally-bound proteins, due to the unique contour geometry 
(Fig. ID). As the adaptive paths wrap around the diagonal contour in the region of high 
Eb and low Ef, the landscape flattens, making selection weaker and substitutions slower 
(Fig. S5C). Thus most of the waiting occurs in the middle of the path rather than the end 
(Fig. S5D). Adaptation accelerates toward the end of the path as the strength of selection 
increases again. If the intermediate slow-down is significant enough, a protein may not have 
time to complete the second half of its path before environmental conditions change, so that 
it will never evolve the new binding function. 

Discussion 

Protein folding and binding as evolutionary spandrels. In the decades since Gould 
and Lewontin's paper [20], the existence of evolutionary spandrels has emerged as a critical 
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evolutionary concept. There are many possible scenarios in which spandrels can evolve [20, 
21], although two key mechanisms are neutral processes, such as genetic drift and biases in 
mutation and recombination [45], and indirect selection arising from coupled traits. Here we 
have focused on the latter, which we expect to be more important on short time scales. 

It has been previously argued that the marginal stability of most proteins may be an 
evolutionary spandrel that evolved due to mutation-selection balance [3, 12, 13]. We sug- 
gest more broadly that having folding stability at all may be a spandrel for proteins with 
no misfolding toxicity. Even more striking is the possibility that some binding interactions 
may be spandrels that evolved solely to stabilize proteins with toxic misfolding; this would 
significantly affect our interpretation of data on proteome-wide interactions [46] . In particu- 
lar, we expect more widespread nonfunctional interactions among proteins with less intrinsic 
stability. Indeed, protein abundance is believed to correlate positively with stability {—Ef) 
to explain the observed negative correlation of abundance with evolutionary rate [18, 19]. 
Furthermore, models of protein-protein interaction networks imply that protein abundance 
also correlates negatively with the number of interactions [16]. Together these argue that sta- 
bility should indeed be negatively correlated with the number of interactions. Experiments 
on specific proteins also support this finding: for example, destabilizing mutations in E. coli 
dihydrofolate reductase were found to be compensated at high temperature by protein bind- 
ing, which protected against toxic aggregation [28]. Previously the role of binding-mediated 
stability has been primarily discussed in the context of intrinsically disordered proteins [47] , 
described by the high Ef regime of our model. 

Pareto optimization of proteins. The Pareto front is a useful concept in problems 
of mult i- objective optimization [44]. The Pareto front in our model consists of the protein 
sequences along the low Ef, low Ef, edge of the genotype distribution (see e.g. Fig. 2A). 
Pareto optimization assumes that all states on the front are valid final states for adaptation; 
this in turn implies that fitness has linear dependence on the individual traits. However, 
nonlinear fitness functions with saturation effects will confound this assumption. Our model 
shows how this nonlinearity leads to a small subset of true final states on or even off the 
front. Thus Pareto optimization does not capture a key feature of the underlying biophysics, 
providing only a rough approximation to the true dynamics. 

Epistasis and evolutionary predictability. Our results also shed light on the role 
of epistasis - the correlated effects of mutations at different sites - in protein evolution. 
Epistasis underlies the ruggedness of fitness landscapes [31, 32]. Magnitude epistasis is 
widespread in our model, while sign epistasis only arises in regions where the fitness contours 
are curved (Fig. 1E,F). This picture is qualitatively consistent with studies of empirical fitness 
landscapes [32] and with directed evolution experiments [3] . 

Epistasis determines the predictability of evolution, an issue of paramount importance 
in biology [25, 33, 34]. In most cases considered here, limited sign epistasis gives rise to less 
predictable intermediate pathways (high 5 pa th) but highly predictable final outcomes (low 

^com ) • 

However, there are two major exceptions to this pattern. First, proteins with a binding 
interaction under no direct selection may have multiple local maxima, some with strong and 
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others with weak binding (Fig. 3). Here both intermediate pathways and final states are 
unpredictable - pure chance, in the form of random mutations, drives the population to one 
binding fate or the other. The second exception occurs in proteins with direct selection for 
both binding and folding. Here there is usually a single maximum, but the adaptive paths are 
tightly constrained in energy space (Fig. 4). Thus evolution of proteins with both functional 
binding and deleterious misfolding, which should include a large fraction of natural proteins, 
is highly predictable at the level of energy traits. 

Methods 

Energetics of protein folding and binding. Folding energetics are probed experimen- 
tally and computationally by measuring the changes in Ef resulting from single point muta- 
tions. Since these changes are observed to be universally distributed over many proteins [8], 
we sample entries of ef from a Gaussian distribution with mean 1.25 kcal/mol and standard 
deviation 1.6 kcal/mol. For the reference sequence a re f, e/(z, a l re{ ) = 0 for alii £ {1, . . . , L}, 
such that Ef(a re f) = EJ { . The parameter E™ n is defined as the binding energy of the geno- 
type <Tbb with the lowest E^\ 65(2,0-^) = 0 for all i G {1,. . . ,1/}. Since binding hotspot 
residues typically have a 13 kcal/mol penalty for mutations away from the wild- type amino 
acid [36, 37], we sample the other entries of 65 from an exponential distribution defined in 
the range of (l,oo) kcal/mol, with mean 2 kcal/mol. This distribution is consistent with 
alanine-scanning experiments which probe energetics of amino acids at the binding inter- 
face [48] . We consider L = 6 hotspot residues and a reduced alphabet of k = 5 amino acids 
(grouped into negative, positive, polar, hydrophobic, and other), resulting in 5 6 = 15625 
possible genotypes. Our population genetics model and the algorithm for exact calculation 
of adaptive path statistics are available in SI Methods. 
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Figure 1: Fitness, selection, and epistasis in energy trait space. (A) Phase dia- 
gram of protein structural states. Dashed lines separate structural phases of the protein 
corresponding to plateaus on the fitness landscape; arrows represent the folding transition 
(green), binding transition (red), and the coupled folding-binding transition (blue). Fitness 
landscapes T{Ef^E^) with direct selection (B) for binding only (/ u b = f u i = 0), (C) for 
folding only (/ u b = 1, f u i = 0), and (D) for both binding and folding (/ u b = 0.9, / u f = 0). 
Black contours indicate constant fitness values. The contours are uniformly spaced in energy 
space; fitness differences between adjacent contours are not all equal. Streamlines indicate 
the direction of the selection "force" V log J 7 , with color showing its magnitude (decreasing 
from red to blue). (E) Projection of a genotype distribution and mutational network into 
energy space for L = 2 and a two-letter (k = 2) alphabet. (F) Blue arrows indicate the 
same mutation on different genetic backgrounds. When the fitness contours are straight, 
the mutation is beneficial regardless of the background {a\ or a 2 ). However, with curved 
contours, the same mutation can become deleterious (a 3 —> 03), indicative of sign epistasis. 
Sign epistasis from curved contours can give rise to multiple local fitness maxima (e.g., AA 
and BB in (E)). 
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Figure 2: Properties of adaptation with direct selection for binding only. (A) Global 
distribution of folding and binding energies for all k L = 5 6 genotypes in a single realization 
of the model with a marginally stable protein (£j ef = —3 kcal/mol). The black star indi- 
cates the initial state for adaptation (global maximum on J^), red triangles indicate local 
fitness maxima on F 2) shaded according to their commitment probabilities ^(cr), and the 
blue crosses indicate best-folding and best-binding genotypes. The magenta line connects 
genotypes on the Pareto front, and the black contours indicate constant fitness Ti- (B) The 
region of energy space accessible to adaptive paths, zoomed in from (A). Example paths are 
shown in blue and green; black circles indicate intermediate states along paths, sized propor- 
tional to their path density ty{p)\ small gray circles are genotypes inaccessible to adaptation. 
(C) Average number m of local fitness maxima (solid, green) and average number m acc of 
local maxima accessible to adaptation (dashed, blue) versus EJ { . The average number of 
maxima is greatest at EJ { w —3 kcal/mol, where multiple local maxima are separated by 
« 2.23 substitutions on average. (D) Average per-residue Hamming distance between the 
maxima and the best-folding genotype (5f] solid, green) and the best-binding genotype 
dashed, blue) versus Ej ef . (E) Average distributions p(t) of path lengths (number of substi- 
tutions) I for stable proteins (EJ { = —15 kcal/mol), marginally stable proteins (Ey ef = —3 
kcal/mol), and intrinsically unstable proteins (Ej ei = 5 kcal/mol). (F) Per-substitution path 
entropy S^^/I (solid, green) and entropy of commitment probabilities S coin (dashed, blue) 
versus EJ { . Panel (E) is averaged over 10 5 realizations of the model; all other averages are 
taken over 10 4 realizations. In all panels / u b = / u f = 0 and E™ n = E™ n = —5 kcal/mol. 
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Figure 3: Properties of adaptation with direct selection for folding only. (A) The 

average number of local maxima m (solid, green) and their average per-residue Hamming 
distance from the best-folding [5f] dashed, blue) and the best-binding {5^ dotted, red) 
genotypes versus E™ m . (B) Probability that adaptation occurs when the binding target is 
changed (i.e., the initial state is not coincident with any of the final states), as a function 
of E™ n and E™ n . (C,D) Example landscape with divergent binding fates: there are two 
accessible local maxima, one with E\> < 0 (favorable binding, i/j(a) = 0.6) and the other 
with Ef, > 0 (negligible binding, i/j(ct) = 0.4). All symbols are the same as in Fig. 2A,B. 
(E) Average distribution of local maxima, weighted by their commitment probabilities. The 
average commitment entropy for realizations with divergent fates is S com ~ 0.43. In (C)-(E) 
we used E™ n = E™ n = —6.5 kcal/mol. (F) The probability of having divergent fates versus 
£™n — E™ n . Panel (E) is averaged over 10 5 realizations of the model; all other averages are 
taken over 10 4 realizations. In all panels / u b = 1, / u f = 0, and £^ ef = 0 kcal/mol. 
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Figure 4: Properties of adaptation with direct selection for both folding and 
binding. (A, B) Distribution of folding and binding energies in an example landscape for a 
marginally stable and marginally bound protein; all symbols are the same as in Fig. 2A,B. 
(C) Landscape averaged over 10 5 realizations. Distribution of initial states is shown in 
green, intermediate states in blue (weighted by their path densities), and final states in 
red (weighted by their commitment probabilities). In all panels / u b = 0.9, / u f = 0, and 
Ef = E™ in = E™ in = -4 kcal/mol. 
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Supplementary Methods 

Population genetics model. In the monomorphic limit, the population is described by 
a single point in genotype space [39]. The population evolves over time via mutations that 
arise sequentially and either fix or disappear. Each fixation event leads to an amino acid 
substitution in the entire population. The rate of making a substitution from genotype a to 
genotype a 1 is given by [35] 

W(a'\a) = Nu </>(*' \a), (SI) 

where TV is the effective population size, u is the mutation rate, and (j>(a'\a) is the probability 
of a single a' mutant fixing in a population of wild-type a. Typically the fixation probability 
depends only on the relative selection coefficient s = J r (cr / )/J r (a) — 1 between the two 
genotypes, where J 7 (a) is the fitness of genotype a. For example, in the Wright-Fisher 
model, (j)(s) = (1 — e _2s )/(l — e~ 2Ns ), where TV is the effective population size [49]. In the 
strong-selection limit (N\s\ 3> 1), 

, , v f 1 - e~ 2s for s > 0 , . 

*W * { 0 for 8 < o (S2) 

Thus the effective population size TV sets the overall time scale (Nu)~ l of substitutions but 
does not affect fixation probabilities. 

Statistics of adaptive paths. We calculate statistical properties of the adaptive paths 
using a transfer matrix-like algorithm [29, 30]. Let S be the set of all genotypes accessible 
to adaptation, and let Sf be the set of final state genotypes (e.g., local fitness maxima). 
Define W^o^o") as the rate of making a substitution from genotype a to genotype a 1 (e.g., 
given by Eq. SI). The rate matrix defines 9(a) = (J^ nn a > of G W(a f \a))~ l , the mean waiting 
time in genotype a before a substitution occurs, where the sum is over all genotypes a 1 one 
mutation away from a (nearest mutational neighbors, "nn"). The substitution rates also 
determine the probability Q(o J \a) = W(a f \a)9(a) of making the substitution a —> a 7 , given 
that a substitution occurs out of a. 
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For each substitution I and intermediate genotype a, we calculate Pt[o\ the total prob- 
ability of all paths that end at a in t substitutions; 7)(cr), the total average time of all such 
paths; and I^(cr), their total entropy. These quantities obey the following recursion relations: 



P t {a>) = Q{^W)Pt-i(<r), (S3) 

nn a of a' 

T t {a>) = E Q(^k)[T<-i(a) + 0(<7)^-i(<7)], 

nn a of a' 

T e (a') = Yl Q(^W)^i-i(<r)-(logQ(<T , \a))P e - 1 (a)], 



nn a of a' 



where Po(cr) = 1 if a is the initial state and Po(&) = 0 otherwise, and T 0 (a) = r 0 (V) = 0 for 
all The final states a G Sf are treated as absorbing to ensure that only first-passage 

paths are counted. We use these transfer matrix objects to calculate the path ensemble 
quantities described in the text: 



aeS f £=1 

* = EE T &) = X>w = E r (^ = Eww. 

A 

^ = E E r <( ff )> r W = £ 

£=1 creS f aes 

The sums are calculated up to a path length cutoff A, which we choose such that 1 — 
Ylf=i P(P) < 10 -6 . Note that the calculations for the state-dependent quantities ip(cr) and 
r(a) are simplified in this model (compared to more general cases [29, 30]) since the strong- 
selection dynamics prevents the population from traversing loops in genotype space. The 
time complexity of the algorithm scales as 0{^N A) [29], where 7 is the average connectivity 
and N is the total size of the state space. For genotypic sequences of length L and an 
alphabet of size fc, 7 ~ L{k — 1) and N ~ k L . 

Validity of the additive energy model. Double mutant experiments indicate that 
the additive energy model is a good approximation for residues that are not in direct physical 
contact [38, 50]. For spatially-close residues, the mutational effects are largely "sub-additive" 
(diminishing-returns magnitude epistasis): two (de)stabilizing mutations combined will still 
usually be (de)stabilizing, but less so than the sum of their individual effects [38, 50]. For 
example, Istomin et al. [50] find that while residues separated by more than 6 A are nearly 
additive (correlation R 2 = 0.97 with a slope of 0.88 between the sum of AAG's for two 
single mutants and A AG for the double mutant), spatially-close residues are substantially 
sub-additive (R 2 = 0.84, slope of 0.54). Nonetheless, in regions with straight contours which 
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represent most of our fitness landscapes, sub-additive energies cannot produce sign epistasis; 
substantial deviations from energy sub-additivity are required to create additional local 
maxima or place significant constraints on adaptive paths. Thus it appears that deviations 
from energy additivity will not lead to qualitative changes in our model's predictions. 
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Figure SI: Average selection strength. (A) Average log 10 s (s is the selection coefficient) 
of all accessible beneficial substitutions as a function of EJ f and E™ n = E™ n = E™ n in 
the case of direct selection for binding only (/ u f = / u b = 0). Due to the symmetry of 
this case (Fig. IB), we can neglect differences in E™ n and E™ n without loss of generality. 

(B) Same as (A) but limited to accessible substitutions that end at local fitness maxima. 

(C) Average log 10 s of all accessible beneficial substitutions as a function of E™ n and E™ n 
in the case of selection for folding only (/ u f = 0, / u b = 1, Ej ei = —5 kcal/mol). (D) Same as 
(C) but limited to accessible substitutions that end at local fitness maxima. (E, F) Same as 
(C, D) but for Ej ef = 0 kcal/mol. Simultaneous selection for both binding and folding yields 
qualitatively similar results. All data points are averages over 10 4 landscape realizations. 
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Figure S2: Example landscapes for stable and intrinsically unstable proteins with 
direct selection for binding only. Symbols and randomly generated energy matrices (ej, 
e& 15 and e^) are the same as in Fig. 2A,B. (A, B) Stable protein (£j ef = — 15 kcal/mol). (C, 
D) Intrinsically unstable protein (£j ef = 5 kcal/mol). As in Fig. 2A,B, f u ^ = / u f = 0 and 
E bt n = = ~ 5 kcal/mol. 
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Figure S3: Average landscapes for direct selection for binding only. As in Fig. 4C, 
the distribution of initial states is shown in green, intermediate states in blue (weighted by 
their path densities), and final states in red (weighted by their commitment probabilities). 
(A) Stable proteins (Ej ei = — 15 kcal/mol). (B) Marginally stable proteins [EJ { = —3 
kcal/mol). (C) Intrinsically unstable proteins (Ej ei = 5 kcal/mol). All landscapes are 
averaged over 10 5 realizations. As in Fig. 2A,B, / ub = / uf = 0 and Ef£ in = £^ in = —5 
kcal/mol. 
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Figure S4: Example and average landscapes for direct selection for folding only. 

Symbols in top and middle panels are the same as in Fig. 2A,B, and the color scheme in 
the bottom panels is the same as in Fig. 4C and Fig. S3. (A) Strong binding to both old 
and new targets (E™ n = E™ n = —8 kcal/mol). (B) Weak binding to old target and strong 
binding to new target (E™ n = 0 kcal/mol, E™ n = —8 kcal/mol). (C) Strong binding to old 
target and weak binding to new target (E™ n = —8 kcal/mol, E™ n = 0 kcal/mol). We use 
fub — 1 5 /uf — 0, and Ej ef = 0 kcal/mol in all cases. In the bottom panels, the landscapes 
are averaged over 10 5 realizations. 
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Figure S5: Distribution of adaptation times over intermediate states. (A) The 

same landscape realization as in Fig. 2A,B (selection for binding only on a marginally stable 
protein), but with each intermediate state a sized proportional to r(cr), the average time 
spent in that state. (B) The probability p(£) (solid, green) of taking an adaptive path of 
exactly t substitutions and the average time r(£) (dashed, blue) spent by paths at the £th 
substitution, averaged over 10 5 realizations with / u b = f u f = 0, EJ { = —3 kcal/mol, and 
^mm _ £jrnm _ _^ kcal/mol. (C, D) Same as (A, B), but with the landscape realization 
used in Fig. 4A,B (selection for both binding and folding, / ub = 0.9, / u f = 0, Ej ei = E™ n = 
E™ n = -4 kcal/mol). 
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